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\ Abstract 

We provide a complete picture to the self-gravitating non-relativistic gas at ther- 
^ \ mal equilibrium using Monte Carlo simulations, analytic mean field methods (MF) 

' and low density expansions. The system is shown to possess an infinite volume limit 

^ in the grand canonical (GCE), canonical (CE) and microcanonical (MCE) ensembles 

when (A/", y) oo, keeping NlV^I"^ fixed. We compute the equation of state (we 
do not assume it as is customary), as well as the energy, free energy, entropy, chem- 
ical potential, specific heats, compressibilities and speed of sound; we analyze their 
properties, signs and singularities. All physical quantities turn out to depend on a 
single variable r\ = ^ that is kept fixed in the iV — > 00 and V ^ oc limit. The 
system is in a gaseous phase for r] < rjx and collapses into a dense object for r] > r]T 
in the CE with the pressure becoming large and negative. At r/ ~ r^^^ the isothermal 
compressibility diverges. This gravitational phase transition is associated to the 
Jeans' instability. Our Monte Carlo simulations yield r]T — 1.515. PV/[NT] = f[rj) 
and all physical magnitudes exhibit a square root branch point at rj = r]c > ?/t- 
The values of rjx and rjc change by a few percent with the geometry for large N: 
for spherical symmetry and N = 00 (MF), we find r]c = 1.561764... while the 
Monte Carlo simulations for cubic geometry yields rjc — 1.540. In mean field and 
spherical symmetry cy diverges as [rjc — rf) ' for rj ] r]c while cp and kt diverge 
as {tjq — r])~^ for r] ] r]Q = 1.51024. . .. The function /(??) has a second Riemann 
sheet which is only physically realized in the MCE. In the MCE, the collapse phase 
transition takes place in this second sheet near rjMC = 1-26 and the pressure and 
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temperature are larger in the collapsed phase than in the gaseous phase. Both 
collapse phase transitions (in the CE and in the MCE) are of zeroth order since 
the Gibbs free energy has a jump at the transitions. The MF equation of state 
in a sphere, /(r/), obeys a first order non- linear differential equation of first kind 
Abel's type. The MF gives an extremely accurate picture in agreement with the 
MC simulations both in the CE and MCE. Since we perform the MC simulations 
on a cubic geometry they describe an isothermal cube while the MF calculations 
describe an isothermal sphere. The local properties of the gas, scaling behaviour of 
the particle distribution and its fractal (Haussdorf) dimension are investigated in 
the companion paper 
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1 Statistical Mechanics of the Self-Gravitating Gas 



Physical systems at thermal equilibrium are usually homogeneous. This is the case for 
gases with short range intermolecular forces (and in absence of external fields). In such 
cases the entropy is maximum when the system homogenizes. 

When long range interactions as the gravitational force are present, even the ground 
state is inhomogeneous. In this case, each element of the substance is acted on by very 
strong forces due to distant particles of the gas. Hence, regions near to and far from the 
boundary of the volume occupied by the gas will be in very different conditions, and, as 
a result, the homogeneity of the gas is destroyed [Q. The state of maximal entropy for 
gravitational systems is inhomogeneous. This basic inhomogeneity suggested us that 
fractal structures can arise in a self-interacting gravitational gas[^ H, H |], 0. 

The inhomogeneous character of the ground state for gravitational systems explains 
why the universe is not going towards a 'thermal death'. A 'thermal death' would mean 
that the universe evolves towards more and more homogeneity. This can only happen 
if the entropy is maximal for an homogeneous state. Instead, it is the opposite what 
happens, structures are formed in the universe through the action of the gravitational 
forces as time evolves. 

Usual theorems in statistical mechanics break down for inhomogeneous ground states. 
For example, the specific heat may be negative in the microcanonical ensemble (not in 
the canonical ensemble where it is always positive) . 

As is known, the thermodynamic limit for self-gravitating systems does not exist in 
its usual form {N —> oo, V —>■ oc, N/V = fixed). The system collapses into a very dense 
phase which is determined by the short distance (non-gravitational) forces between the 
particles. 

We instead find that the thermodynamic functions exist in the dilute limit 

A^ oo , V ^ oo , yYj^ = fixed (1) 

where V stands for the volume of the box containing the gas. In such a limit, the energy 
E, the free energy and the entropy turns to be extensive. That is, we find that they take 
the form of A^ times a function of 

Gm'^N , EL 

r] = or 4 



LT ^ Gm?N'^ 

where 'q and ^ are intensive variables. Namely, rj and ^ stay finite when A^ and V = Lr" 
tend to infinite, r] is appropriate for the canonical ensemble and ^ for the microcanonical 
ensemble. Physical magnitudes as the specific heat, speed of sound, chemical potential 
and compressibility only depend on 77 or ^. rj and ^ as well as the ratio N/ L are therefore 
intensive magnitudes. The energy, the free energy, the Gibbs free energy and the entropy 
are of the form A^ times a function of rj. These functions of rj have a finite N = 00 limit 
for fixed rj (once the ideal gas contributions are subtracted). Moreover, the dependence 
on 1] in all these magnitudes express through a single universal function /(//). 

We study here and in the companion paper (called paper II in what follows) the sta- 
tistical mechanics of the self-gravitating gas. That is, our starting point is the partition 
function for non-relativistic particles interacting through their gravitational attraction 
in thermal equilibrium. We study the self-gravitating gas in the three ensembles: micro- 
canonical (MCE), canonical (CE) and grand canonical (GCE). We performed calculations 
by three methods: 
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• By expanding the partition function through direct calculation in powers of 1/^ 
and 1] for the MCE and CE, respectively. These expressions apply in the dilute 
regime {C, ^ 1 , r] ^ 1) and become identical for both ensembles for iV — > oo. At 
77 = 0=1/^ we recover the ideal gas behaviour. 

• By performing Monte Carlo simulations both in the MCE and in the CE. We found 
in this way that the self-gravitating gas collapses at a critical point which depends 
on the ensemble considered. As shown in fig. 1 the collapse occurs first in the 
canonical ensemble (point T) . The microcanonical ensemble exhibits a larger region 
of stability that ends at the point MC (fig. 1). Notice that the physical magnitudes 
arc identical in the common region of validity of both ensembles within the statistical 
error. Beyond the critical point T the system becomes suddenly extremely compact 
with a large negative pressure in the CE. Beyond the point MC in the MCE the 
pressure and the temperature increase suddenly and the gas collapses. The phase 
transitions at T and at MC are of zeroth order since the Gibbs free energy has 
discontinuities in both cases. 

• By using the mean field approach we evaluate the partition function for large A'^. We 
do this computation in the grand canonical, canonical and microcanonical ensem- 
bles. In the three cases the partition function is expressed as a functional integral 
over a statistical weight which depends on the (continuous) particle density. These 
statistical weights are of the form of the exponential of an 'effective action' propor- 
tional to A^. Therefore, the A" — > 00 limit follows by the saddle point method. The 
saddle point is a space dependent mean field showing the inhomogencous character 
of the ground state. Corrections to the mean field are of the order 1/A^ and can be 
safely ignored for N ^ 1 except near the critical points. These mean field results 
turned out to be in excellent agreement with the Monte Carlo results and with the 
low density expansion. 

We calculate the saddle point (mean field) for spherical symmetry and we obtain from 
it the various physical magnitudes (pressure, energy, entropy, free energy, specific heats, 
compressibilities, speed of sound and particle density). Furthermore, we compute the 
determinants of small fluctuations around the saddle point solution for spherical symmetry 
for the three statistical ensembles in paper 11. 

When any small fluctuation around the saddle point decreases the statistical weight 
in the functional integral, the saddle point is dominating the integral and the mean fleld 
approach is fully valid. In that case the determinant of small fluctuations is positive. A 
negative determinant of small fluctuations indicates that some fluctuations around the 
saddle point arc increasing the statistical weight in the functional integral and hence the 
saddle point does not dominate the partition function. The mean field approach cannot 
be used when the determinant of small fiuctuations is negative. 

The zeroes of the small fluctuations determinant determine the position of the critical 
points for the three statistical ensembles. The Monte Carlo simulations for the CE and 
the MCE show that the self-gravitating gas collapses near the critical points obtained 
from mean field. 

The saddle point solution is identical for the three statistical ensembles. This is not 
the case for the fluctuations around it. The presence of constraints in the CE (on the 
number of particles) and in the MCE (on the energy and the number of particles) changes 
the functional integral over the quadratic fluctuations with respect to the GCE. 
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The saddle point of the partition function turns out to coincide with the hydrostatic 
treatment of the self-gravitating gas|^-||TB[. (Usually known as the 'isothermal sphere' in 
the spherically symmetric case). 

Our Monte Carlo simulations are performed in a cubic geometry. The equilibrium 
configurations obtained in this manner can thus be called the 'isothermal cube'. 

We find for spherical symmetry: tjqq = 0.797375 . . . , r]^ = 2.517551 . . . and rj^jc ~ 
2.03085 . . .. The variable t]^ appropriate for a spherical symmetry is defined as t]^ = 
= ^ (^)^^^ = 1.61199 ... r]. 

The values of rjT and t]c change by a few percent with the geometry and with the 
number of particles (for large > 500). For spherical symmetry and N = oo (mean 

field) we obtain rjc = ( J;:) ^ Vc ~ 1-56176 . . .. Our Monte Carlo simulations yield 
rjT — 1.515. We find from the mean field approach that the isothermal compressibility 
diverges ai r] = rj^ = 1.51024 . . . c:^ tjt for spherical symmetry. 

The conclusion being that the MF correctly describes the thermodynamic limit except 
near the critical points (where the small fiuctuations determinant vanishes); the MF is 
valid for N\r] — rjcritl S> 1. The vicinity of the critical point should be studied in a double 
scaling limit ^ oo, r] ^ rjcrit- 

In summary, the picture we get from our calculations using these three methods show 
that the self-gravitating gas behaves as a perfect gas for 77 — >• 0, 1/^ — > 0. When 77 and 
1/^ grow, the gas becomes denser till it suddenly condenses into a high density object at a 
critical point GC, C or MC depending upon the statistical ensemble chosen. In the Monte 
Carlo simulations for the CE the collapse takes place at the point T slightly before rjc. rj 
is related with the Jeans' length dj of the gas through 77 = 3(L/(ij)^. Hence, when r] goes 



beyond ?7t, the length of the system becomes larger than dj/^Jriq-fS. The collapse at T in 
the CE is therefore a manifestation of the Jeans' instability. The saddle point ceases to 
describe the physics at C since the determinant of fiuctuations for the CE vanishes there. 

In the MCE, the determinant of fiuctuations vanishes at the point MC. The physical 
states beyond MC are collapsed configurations shown by the Monte Carlo simulations 
[see fig. Actually, the gas collapses in the Monte Carlo simulations slightly before the 
mean field prediction for the point MC. The phase transition at the critical point MC is 
the so called gravothermal catastrophe fT^ . 



The gravitational interaction being attractive without lower bound, a short distance 
cut-off {A) must be introduced in order to give a meaning to the partition function. We 
take the gravitational force between particles as —G m^/r^ for r > A and zero for r < A 
where r is the distance between the two particles. We show that the cut-off effects are 
negligible in the N = 00 limit. That is, once we set N = 00 with fixed t], all physical 
quantities are finite in the zero cut-off limit {A = 0). The cut-off effects are of the order 
A'^/L'^ and can be safely ignored. 

All physical quantities are expressed in terms of f{r]). Besides computing f{ri) nu- 
merically in the mean field approach, we obtain analytic results about it from the Abel's 
equation. There is a square root branch point in /(?]) at rjc. The specific heat is positive 
in the first sheet and negative in the second sheet. This second sheet is only physically 
realized in the microcanonical ensemble (MCE). [The specific heat is positive definite 
in the CE]. f{ri) has infinitely many branches in the rj plane but only the first two are 
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physically realized. Beyond MC the states described by the mean field saddle point are 
unstable. We plot and analyze the equation of state, the energy, the entropy, the free 
energy, cy, cp, the isothermal compressibility and the speed of sound [figs. p!0| - |T5[| . Most 



of these physical magnitudes were not previously computed in the literature as functions 

of T]. 

We find analytically the behaviour of f{ri) near rjc in mean field, 

fMPiv"") " I + 0.213738 . . . ^vE-v"" + 0.172225 . . . (ry^ - r/^) + O [(r^^ - r^y/' 

This shows that the specific heat at constant volume diverges as {rjc — v)^^^^ T Vc- 

The specific heat at constant pressure and the isothermal compressibility diverge at rjo as 
(^0 ~ v)~^- These mean field results apply for \ri — rjcl ^ 1 ^ N\ri — r]c\. Fluctuations 
around mean field can be neglected in such a regime. 

The Monte Carlo calculations permit us to obtain f{ri) in the collapsed phase. Such 
result (which is cutoff dependent) cannot be obtained in the mean field approach. The 
mean field only provides information (as f{f])) in the gas phase. 

For the self-gravitating gas, we find that the Gibbs free energy $ is not equal to 
times the chemical potential and that the thermodynamic potential Q is not equal to 
—PV as usual|0]. This is a consequence of the dilute thermodynamic limit N oo, L 
oo, N/L =fixed. 

We compute local properties of the gas in paper II. That is, the local energy den- 
sity e(r), local particle density and local pressure. Furthermore, we analyze the scaling 
behaviour of the particle distribution and its fractal (Haussdorf) dimension. 

This paper is organized as follows. In section II we present the statistical mechanics 
of the self-gravitating gas in the microcanonical ensemble, in sec. Ill we do the analogous 
presentation for the canonical ensemble and in sec. IV we contrast the results for the CE 
and the MCE. Sec. V contains the results from Monte Carlo simulations and we develop 
in sec. VI the mean field approach. In sec. VII we present the results for intensive 
magnitudes. Discussion and remarks are presented in section VIII whereas appendixes 
A-C contain relevant mathematical developments. 

2 Statistical Mechanics of the Self-Gravitating Gas: 
the microcanonical ensemble 

We investigate in this section an isolated set of non-relativistic particles with mass m 
interacting through Newtonian gravity with total energy E. That is, a self-gravitating 
gas in the microcanonical ensemble. We assume the system being on a cubic box of side 
L just for simplicity. We consider spherical symmetry in sec. VI. Please notice that we 
never use periodic boundary conditions. 

At short distances, the particle interaction for the self-gravitating gas in physical situ- 
ations is not gravitational. Its exact nature depends on the problem under consideration 
(opacity limit. Van der Waals forces for molecules etc.). We shall just assume a repulsive 
short distance potential, that is, 

1 f ~wki\ for > A 

VAm-qj\) = - \^ = (2) 

l^^-^^l^ [ +i ioT\q,-q,\<A 
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where A « L is the short distance cut-off. 

The presence of the repulsive short-distance interaction prevents the collapse (here 
unphysical) of the self-gravitating gas. In the situations we are interested to describe 
(interstellar medium, galaxy distributions) the collapse situation is unphysical. 

The entropy of the system can be written as 



S{E,N)^ log 



1 1 / /A^' 



6 



N 2 



(3) 



where 



i<i<j<N\1i ^iU 



(4) 



and G is Newton's gravitational constant. 

In order to compute the integrals over the momenta pi, (1 < / < N), we introduce 
the variables, 



1 



'2rn 



Pi 



We can now integrate over the angles in 3A^ dimensions. 



N J3 



r+oo r+oo - g-^p^ 

y_oo ■■■i-oo fJi (27r)3 

iN/2 ^oo 

IZT Jo ' 



N 



E-J2Pi-U{qu...qN) 

1=1 



^3iV-l 



27r 



m 



\ 3N/2 



r(f) 



dps E- -U{qi, . . .^jv) 



27ry 



r(f) 



[E - U{q^, . . . qN)]'"^'-' e [E - U{q„ . . . q^)] 



(5) 



The delta function in the energy thus becomes the constraint of a positive kinetic energy 
E — U{qi, . . . qjy) > 0. We then get for the entropy, 

{( ^^^"^ L ^ 

;;;V3^X / •••/ Y[d\[E-U{q,,...q^)r"-'e[E-U{q,,...q^)] 
\T ) '=1 

(6) 

It is convenient to introduce the dimensionless variables fi, 1 < I < N making explicit 
the volume dependence as 



qi ^ Lri , fi = {xi,yi,zi) , 
< xi,yi,zi<l. 



(7) 



That is. in the new coordinates the gas is inside a cube of unit volume. 
The entropy then becomes 



S{E,N) = log 



J\f3N-2 ^97V/2-2 ^37V/2+l Q3N/2-1 



7V!r(f) {2n) 



3N/2 



d\ 



3N/2-1 



(8) 
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where we introduced the dimensionless variable ^, 

EL 



and 



rN) = 



Gm? m 
1 



E 



1 



ri - r,-L 



(9) 
(10) 



l<l<j<N 

where a = A/ L 1. 

Let us define the coordinate partition function in the microcanonical ensemble as 



Therefore, 



N 



n d'ri 



1=1 



^ + — M(ri,...,rjvJ 



3iV/2-l 
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^ + — M(fi,...,f7v) 



'111 



SiE,N) = log 



J\f3N-2 j^3N/2+l Q3N/2-1 



m 

2 



(2vr) 



3N/2 



+ log w{^,N) . 



We can now compute the thermodynamic quantities, temperature and pressure through 
the standard thermodynamic relations 



1 
T 



'dS\ 
dEj 



and p = T 



V 



dV 



;i2) 



where V = stands for the volume of the system and p is the external pressure on the 
system. 

We obtain the temperature as a function of E and ^ from eqs.(p|) and (|12|) 



T Ed^ ^ ' 



3iVe 
2E 



2 

3N. 



< 



1 



> 



(13) 



where 



< 



u[ri, ...,rNj 



3N/2-2 



,rNj 



So ■■■Si Till d^Ti [e + i«(fi,...,r-W) 



3N/2-1 



,rN, 



The equation of state follows from eqs.(H) and ( [121) 



NT 2 3N 3Nd^ ^ ' 2 



1 + 



3Nj 2 e + 



> 



2 

3iV, 



(14) 



(15) 



We are interested in the large size limit where — > oo, L ^ oo and E oo. We 
consider that ^ = ^ ^2^jv^ stays fixed in such limit. That is, we assume E/N and L/N 
bounded and nonzero when E, L and N —>■ oo. We shall see below that such limit is 
meaningful. 

It is possible to write the energy and the equation of state in terms of a single function 



> 



1 



2 

3N 



(16) 
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We find from eqs.(|g), (0) and 

pV 

iVT 
E 



NT 



1 1 

2 ^ 3^^^'' ^ 3iV 



(17) 



We obtain the virial theorem by ehminating g{^) in the eqs.(|T^ 

pV 



NT ^E + T 



In the case of a perfect gas (no gravity) we have u{.) = 0, g{^) = |, pV = NT and 
= |A^T as it must be. 



where the term T/3 can be neglected for large A^. 
E 

The function g[^) is computed by Monte Carlo simulations, mean field methods and, 
in the weak field limit ^ << 1, is calculated analytically in powers of 1/.^ in subsection 
II.A. 

The specific heat per particle is given by 



cv 



T /dS\ 
N [dfj 



V 



Hence, using eq.(0) yields 



Cv 



5(C) 



or 



1 

Cv 



d_ 



We can relate the specific heat cy with the fiuctuations as follows. We can express 

i + ^^(fi. 



g{^) as an average value using eqs.(|T4D and ([T6| ) 



ZN/2-l 



9 



9i^) ^-nS',...SIY{Ii d^Ti [e+^w(ri,...,f^ 
Computing the derivative with respect to ^ yields. 



ZN/2-2 
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^ + ^M(fi, 



1 

Cv 



2 
3 



A 



< 



where 



A- 



= N{< 



1 



> - < 



(19) 



> 



is of order N^ for ^ oo. [Notice that in the calculation of the fiuctuations we must 
keep the 1/A^ corrections till the end]. 

We can express cy in terms of the fiuctuations of the inverse temperature P = 1/T 
using eq.(|l^): 



1 

Cv 



2 

3 



(20) 



It must be noticed that in the microcanonical ensemble, cy may be positive as well as 
negative. In fact, it becomes negative when the fluctuations are large enough [see sec. V 
and VI]. 

We see that extensivity holds here in an specific way. T, S/ N and pV/ N are of order 
one for — ^ oo provided ^ stays fixed in such limit. That is, we must keep E/N and 
L/N fixed in the N ^ oo limit. 



2.1 The diluted regime: ^ >> 1 

We can obtain the thermodynamic quantities as a series in powers of 1/.^ just expanding 
the integrand in eq.(|Tip- 
We find 



w 



+ 



2^ 



^2 V 3NJ V 3NJ 



where bo, hi and 62 are pure numbers, 

bo -- 
hi - 



+ oirn (21) 



6 Jo ^0 In — 
1 d^ri S'r2 d^r3 
Jo \ri - r2||fi - fal 



1 rl 



(Pri d?r2 



62 = / / 1^ -.2 (22) 
Jo Jo ri — r2 r 



For the cubic geometry chosen, it takes the value 

&r = 1 /V - ^) dx [\l - y) dy f , = 0.31372 , 

o Jo Jo Jo + + -2^ 

For a sphere of unit volume we find 



^sphere ^ _ _ ^ 0.32239839 . . . 



3 / 

^sphere ^ ^ fl^^V^' = 3.786412026 



35 V 3 / 

^s^phere ^ ^( — ) = 5.846665629.... (23) 



4 V 3 / 

We see that the coefficient ho for cubic and spherical geometries only differ by about 

3%. 

We thus find from eq.(^) in the N ao limit 

lim^^o^i log^(e, N) = l\og^+^-^ + ^^{hi- ml) + 0{r') (24) 
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We considered here these integrals in the zero cutoff limit since bo, bi and 62 have 
finite zero cutoff limits. It is easy to see that their finite cutoff values behave as 

boia) - bo = 0{a^) , 61(a) -61 = O(a^) , 62(a) - 62 = 0(a) (25) 



Inserting eq.(p4|) into eq.(|T6|) yields 

3 9 60 



9{0 = l-'-^-^,{h-42bl) + O{n (26) 



and 

pV ^ 3bo 3 



1 - 



{b,-A2bl)+OiC') ■ (27) 



NT ' 2^ 4^2 

We see that after letting ^ 00 the zero cutoff limit is finite. We further discuss this 
important issue in the next section. 

3 Statistical Mechanics of the Self-Gravitating Gas: 
the canonical ensemble 

We investigate in this section the self-gravitating gas considered in the previous section 
but in thermal equilibrium at temperature T = That is, we work in the canonical 
ensemble where the system of N particles is not isolated but in contact with a thermal 
bath at temperature T. We keep assuming the gas being on a cubic box of side L. 
The partition function of the system can be written as 



1=1 



(27r)3 



where 

M 

Pi r^™2 



Zc{N, T) = ^ / . . . / n e-^^- (28) 

N 2 

HN = T.li--Grn^ E (29) 

1=1 i</<i<A^ m-ijiA 



G is Newton's gravitational constant. 

Computing the integrals over the momenta pi, (1 < / < N) 



d^p _ft 
e 2 



(27r)3 \2'k(3) 



yields 

— TV 



We make now explicit the volume dependence introducing the variables r;, 1 < / < 
defined in eq.(|^). The partition function takes then the form, 



3iV 



Zc{N,T) = 4 ' f\.. f'Y[d\ e^'^^-'-'-) , (31) 



A^! V 27r J Jo Jo 



1=1 
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where we introduced the dimensionless variable rj 

Gm'^N 



V 



LT 



and u{fi, . . . , r/v) is defined by eq. (|TO|) . Recall that 

Gm^N 



U 



L 



U{fi, . . . ,fAr) 



(32) 



(33) 



is the potential energy of the gas. 

The free energy takes then the form, 



-T log ZciN,T) = -NT log 



where 



<l>7v(?7) = log 



1 N 



^ri u{ri,...,fN) 



(34) 



(35) 



1=1 



The derivative of the function ^n{v) ^ill be computed by Monte Carlo simulations, mean 
field methods and, in the weak field limit rj << 1, it will be calculated analytically. 
We get for the pressure of the gas, 

[dV , 



P 



NT vT 



(36) 



[Here, V = L"^ stands for the volume of the box and p is the external pressure on the 
system.]. We see from eq. (|35|) that ^n{v) increases with t] since u{.) is positive. Therefore, 
the second term in eq.(^) is a negative correction to the perfect gas pressure 
The mean value of the potential energy U can be written from eq.( 

<U>= -Tr] $'^(7/) 

Combining eqs. (|36[) and ( ^7|) yields the virial theorem. 



as 



pV 
NT 



^ <U > pV 
1 H — „ or 



3NT 



NT 



1 E 

2 ^ 3NT 



(37) 



(3^ 



where we use that the average value of the kinetic energy of the gas is ^NT. 
A more explicit form of the equation of state is 



NT 3N 



(39) 



where 



-<S>n(v) 



r;u{ri,...,riv) 



-(A^-1) e 



/ . . . / T\d^ri M(fi, . . . ,fAr) e 
Jo Jo iJ[ 

Jo Jo f-J^ \ri -rsia 



-$jv(»?) 



^r;u(ri,...,rjv) 



(40) 



This formula indicates that ^'^ir]) is of order for large A^. Monte Carlo simulations as 
well as analytic calculations for small 1] show that this is indeed the case. In conclusion, 
we can write the equation of state of the self-gravitating g 

pV 



NT 



fiv) , 



(41) 
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where the function 



fiv) = 1 - ^ ^Niv) 



is independent of for large and fixed rj. [In practice, Monte Carlo simulations show 
that f{ri) is independent of for > 100]. 
We get in addition, 

<U>= -3NT [1 - firj)] . (42) 
In the dilute limit, ?7 — and we find the perfect gas value 



/(O) = 1 . 



Equating eqs. (|39|) and ( PD yields. 



X 



Relevant thermodynamic quantities can be expressed in terms of the function /(//). We 
find for the free energy from eq. (|34D , 



F = F.-3NT / dx 



V l-f{x) 



X 



where 



Fo = -NT log 



eV /mT\3/2' 
aT \ 2^) 



(43) 



(44) 



is the free energy for an ideal gas. 

We find for the total energy E, chemical potential /i and entropy S the following 
expressions. 



E = 3NT 



fiv) 



(45) 



fx 

S 



(OF) 



-3T[l-f{r])]-3T f'dx 

Jo 



X 



V 



So + 3N 



'dx^-^^ + fiv)-l 



X 



(46) 



where 



^n = -5 + -Ar. 



T 2 



is the entropy of the ideal gas. 

Notice that here the Gibbs free energy 

^ = F + pV = Fq + NT 



X 



(47) 



is not proportional to the chemical potential. That is, here ^ fiN and we have instead, 

^ - fiN = 2NT[1- f{r])] . (48) 
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This relationship differs from the customary one (see p[) due to the fact that the dilute 
scaling relation ~ L holds here instead of the usual one ~ L^. The usual relationship 
is only recovered in the ideal gas limit = 0. 



The specific heat at constant volume takes the forml 



N [df. 



fiv) - V fiv) 



(49) 



where we used eq.(^). This quantity is also related to the fluctuations of the potential 
energy (Af/)^ and it is positive defined in the canonical ensemble, 



cy = - + {Auy 



Here, 



(AUf 



iVT2 



3 [fiv) - V fiv) - 1] ■ 



The specific heat at constant pressure is given by 

fdp\ ■ 

\dv)T 



Cp = Cy 



(50) 
(51) 

(52) 



and then, 



Cp 



Cy + 



[fiv)-vf'iv)f 
fiv) + Ivf'iv) 
^fiv) [fiv)-vf' 



_3 

The isothermal (-R't) and adiabatic iKs) compressibilities take the form 

I =JL I 

\ iv r /(,,) + !,,/'(,,) • 

It is then convenient to introduce the compressibilities 

1 



(53) 







(dV 


Kt = 




\ dp 






(dV 


Ks = 


-V 


\ dp 



Kt = — n— Kt 



and Ks = —TT- K s 



Cy 



Kt 



V fiv) + Ivf'iv) V CP 

which are both of order one (intensive) in the N , L ^ oo limit with N/ L fixed. 



(54) 



(55) 



The speed of sound Vg can be written as ||T9| 



Cp V"^ ( dp 



cv N \dV 



N 



T 



dp ' 
Ncv [df 



where we used eq. (^2]) in the last step. Therefore, 



T 



[fiv)-vf'iv)f 



V 



1 



' dp ' 
dV 



(56) 



fiv) - vf'iv) 



+ fiv) + -^vf'iv) ■ 



(57) 
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The pressure p used in this calculation corresponds to the pressure on the surface of the 
system. Hence, this is the speed of sound on the surface of the system, this is different 
from the speed of sound inside the volume since the ground state is inhomogeneous. We 
compute the speed of sound as a function of the point in paper II. 

We see that the large N limit of the self-gravitating gas is special. Energy, free energy 
and entropy are extensive magnitudes in the sense that they are proportional to the 
number of particles N (for fixed rj). They all depend on the variable r] = '^^^J^ which 
is to be kept fixed for the thermodynamic limit (A^ oo and V oo) to exist. Notice 
that 7] contains the ratio N/L = N V~^^^ which must be considered here an intensive 
variable. Here, the presence of long-range gravitational situations calls for this new 
intensive variable in the thermodynamic limit. 

In addition, all physical magnitudes can be expressed in terms of a single function of 
one variable: f{ri). 



3.1 The diluted regime: 77 << 1 



We can obtain the thermodynamic quantities as a series in powers of rj just expanding 
the exponent in the integrand of ^n{v) [eq-(|35D]. 
To first order in rj we get. 



$jv(r/) 



r]{N -1) 



1 N 

1=1 

1 '■1 d^ri d^r2 



u(ri, 



Jo ri - a 



= 3{N -l)bo7] + 0{7] a^) + 0{Tf) . 

where the coefficient 60 is defined by eq . (^) . 

To first order in rj we see that the cutoff effect is negligible ~ O{o?) [see (p5|)]. 

To second order in rj we find from eq. (|35|) , 



(5^ 



3*Jv(r?) 



N 



\{d\ 



^rj n(ri,...,rjv) 



1=1 

l + 3{N-l)boV 



+ 



2 m 



N{N -1){N- 2){N -3) f d^n d^r^ d^r^ 



\ri - r2 rs - r4 



+ A^(iV-l)(A^-2) 



d^Ti d^r2 d^r^ 



+ 



A^(A^ 



d^ri d^r2 

Ifi - 



(59) 



where the coefficients in front of the integrals count the number of combinations of par- 
ticles yielding the same contribution. Using the notation defined by eqs.(p3) we get 



l + 3(iV- 1) hor] + 7f 



9(iV-i)(iV-2)(iV-3) 2 
2N 
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(iV-2)(iV-2) (iV-2) ■ 

2N ^ 4Ar 2 

Taking the log we get in the infinite hmit; 



(60) 



+ 0{rf) . 



where we have now set a = 0. 

The cutoff effect is here again of order ~ O{o^). It must be noticed that the coefficient 
1)2 which has the stronger dependence on the cutoff [see (|25|) ] cancels out in the N = oo 
limit. 

We therefore find in the low density and the large N limit using eqs.(B^), (^Tj) and 



pV 
iVT 



-bi - I2bl 



+ 0{rf 



(61) 



Furthermore, the speed of sound approaches for — its perfect gas value, 



^io 5 



T 



3 3 ' 9 ' 



6i-36 62 +0{ri^ 



where we used eqs. (|57|) and (|6TD. 

As we see, there are no divergent contributions in ^Niv) the zero cutoff limit to the 
second order in rj. 

At order three a logarithmically divergent integral appears in e*'^'^''-'. Namely, 



3! 



\N(N-l)j 



d^Ti d^r2 if 

~ — log a 



|ri - r2|^ 



A^ 



This integral gives to /(r^) and the other physical magnitudes a contribution of the order 

Therefore, such quantities can be safely neglected for N ^ oo and fixed (small) a since 
f{ri) is of order for A^ — * oo. 

More generally, to the nth. order in rj and n > 3 the leading divergent contribution to 
g*Ar(»7) for a — is of the form 



n\ A^'^2 



N{N - 1) 



In - f2\^ 



T]" 



-,3— n 



! A^"-2 



n 



This gives to f{ri) and the other physical magnitudes a contribution of the order 



1]'' 



-,3—n 



! A^"-i 



n 



As in the n = 3 case, such contributions are negligible in the N ^ oo limit since we take 
it at fixed (small) a. 
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4 Microcanonical vs. Canonical Ensembles 



Let us compare the thermodynamical quantities computed in the microcanonical and 
canonical ensembles in the N oo limit keeping ^ and 1] fixed, respectively. 

We consider here the dilute limit where we dispose of analytic expressions. The Monte 
Carlo and mean field results for the two ensembles will be compared in the next sections 
and in paper II. 

In the dilute limit, we have the expressions ( P7D and ( pT| ) for the equation of state 
in the microcanonical and canonical ensembles, respectively. We want to know whether 
they are or not equivalent. 

Let us start from the microcanonical equation of state (0). We have to express rj in 
terms of C, in order to compare with the canonical equation of state (|6T|) . 

It follows from eqs.®, and (0) that 



Hence, for large ^ and small rj 

V 

and then. 



3 9 bo 



1 



V 



(62) 



(63) 



One easily sees that inserting eq.(|63D in the microcanonical equation of state (p7|) yields 
the canonical equation of state (^Tj) [up to 0{rf) = (9(^~^)]. 

Conversely, starting from the canonical ensemble, it follows from eqs.([l] 
(0) that 



E 

NT 



9(0 



m - 2 



--3boV-V^ [bi-3Qbl]+0{r]^) 



(^) and 
(64) 



and 



V 



m - 2 



3 

2t] 



l-2bo'n--rf (6i - 36 6g) + (9(7/3) 



We see that this relation is identical to eqs. (p^) and ( |53D obtained in the microcanonical 
ensemble [up to 0{rf) = 0{^^^)]. 

Inserting now eq.(|5^ into the canonical equation of state (El) yields the microcanon- 
ical equation of state (|3) [up to 0{r]^) = 0{^-^)]. 

One verifies in the same way that all thermodynamical quantities coincide to the same 
order in both ensembles. 

In summary, the microcanonical and canonical ensembles yield the same results for 

^ cxD to the orders rj'^,-)] and r/^ (or equivalently ^^,C,~^ and ^^'^). 



5 Monte Carlo calculations 

We have applied first the standard Metropolis algorithm [pO[| to the self-gravitating gas 
in a cube of size L in the canonical ensemble at temperature T. We computed in this 
way the pressure, the energy, the average density, the potential energy fiuctuations, the 
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average particle distance and the average squared particle distance as functions of rj. 
We implement the Metropolis algorithm changing at random the position of one particle 
chosen at random. The energy of the configurations is calculated performing the exact 
sums as in eg. ([TI1|) . We used as statistical weight for the Metropolis algorithm in the 
canonical ensemble, 

77 n(ri,...,rjv) 

1 

which appears in the coordinate partition function (|35|) . 

The number of particles went up to 2000. We introduced a small short distance cut- 
off A = 10~^L — 10~^L in the attractive Newton's potential according to eq.(0). All results 
in the gaseous phase were insensitive to the cutoff value. The partition function calcula- 
tion turns to be much less sensible to the short distance singularities of the gravitational 
force than Newton's equations of motion for particles. That is, solving the classical 
dynamics for A^ particles interacting through gravitational forces as well as solving the 
Boltzman equation including the A^-body gravitational interaction requires sophisticated 
algorithms to avoid excessively long computer times [O. As is clear, solving the A^-body 



classical evolution or the kinetic equations provides the time-dependent dynamics and out 
of thermal equilibrium effects which are out of the scope of our approach. 

In the CE, two different phases show up: for rj < rjjn we have a non-perfect gas and 
for ?7 > 772- it is a condensed system with negative pressure. The transition between the 
two phases is very sharp. This phase transition is associated with the Jeans instability. 

A negative pressure indicates that the free energy grows for increasing volume at 
constant temperature [see eq.(p6D]. Therefore, the system wants to contract sucking on 
the walls. 

We plot in figs. 1 and 2, /(r^) = pV/[NT] and (Af/)^ as functions of rj, respectively. 

We find that for small t], the Monte Carlo results for pV/[NT] well reproduce the 
analytical formula (|6TD. pV/[NT] monotonically decreases with rj. 

In the Monte Carlo simulations the phase transition to the condensed phase happens 
for rj = rjT slightly below rjc- For A^ = 2000 we find riT ~ 1.515. For rjT < rj < rjc, the 
gaseous phase can only exist ClS Sb metastable state. 

The average distance between particles < r > and the average squared distance 
between particles < > monotonically decrease with rj. When the gas collapses at 
T]T ,< r > and < > exhibit a sharp decrease. 

The values of pV/[NT], < r > and < > in the condensed phase are independent 
of the cutoff for a < 10~^. The Monte Carlo results in this condensed phase can be 
approximated for 77 > 2 as 

vV 

j— = f(r])^l-K7] , <r>^ 0.016. (65) 

where K ~ 14. 

Since /(r/) has a jump at the transition, the Gibbs free energy $ is discontinuous and 
we have a phase transition of the zeroth order. We find from eq. (^Tj) 

t(colla.p^e)^- i(r,r) ^ ^^^^^^^^^^^ _ ^^^^^ ^ ^ ^ 
We can easily compute the latent heat of the transition per particle (g) using the fact 



that the volume V stays constant. Hence, q = AE/N and we obtain from eq.(^5D 

I = ^'^"""P;^' = 3 [/(collapse) - fM] ^2-3K ^ -62 < . (67) 
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Figure 1: f{f]^) = PV/[NT] as a function of t]^ by Monte Carlo simulations for the 
microcanonical and canonical ensembles (A^ = 2000). Both curves coincide within the 
statistical error till the point T. 
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3 - 



2 - 



1.41 




1.51 



Figure 2: (Af/)^ = 



<u ^>-<u >^ 

N T2 



3 [fiv) ~ V f'iv) — 1] as a function of r] in the gaseous 



phase from Monte Carlo simulations with 2000 particles in the canonical ensemble. Recall 
that cv = 3/2 + {AUf. 



This phase transition is different from the usual phase transitions since the two phases 
cannot coexist in equilibrium as their pressures are different. 

Eq. (|65|) can be understood from the general treatment in sec. Ill as follows. We have 
from eqs.(|D-(|OD 

fiv) = l-l<-> . (68) 
6 r 

The Monte Carlo results indicate that < ^ >~ 42 is approximately constant in the 
collapsed region as i 
such value of < - >. 



collapsed region as well as < r > and < >. Eq.(|03D thus follows from eq. (pHD using 



The behaviour of pV/[NT] near rjc in the gaseous phase can be well reproduced by 



f{ri) = fc + A V^J^ (69) 

where fc ^ 0.316, A ~ 0.414 and r]c ~ 1.540. 

In addition, the behaviour of (Af/)^ in the same region is well reproduced by 

(AUV ^tgo CA- (70) 

Vvc - V 

with C ~ -1.64 and D ^ 0.901. [Notice that for finite N, {AUy will be finite albeit very 
large at the phase transition]. Eq.(^Il) relating /(r/) and (Af/)^ is satisfied with reasonable 
approximation. 
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We thus find a critical region just below rjc where the energy fluctuations tend to 
inflnity as rj ^ rjc- 

The point riT where the phase transition actually takes place in the Monte Carlo 
simulations is at 77^ — 1-51 < rjc- This value for tjt is close to the point rj^ where the 
isothermal compressibility diverges (see sec. VII). They are probably the same point. 

Since Monte Carlo simulations are like real experiments, we conclude that the gaseous 
phase extends from rj = till rj = 7]t i^ti the CE and not till r] = rjc- Notice that in the 
literature based on the hydrostatic description of the self-gravitating gas |l^, |15|, |T^, 



only the instabihty at rj = rjc is discussed whereas the singularities at rj = rjo are not 
considered. 

We then performed Monte Carlo calculations in the microcanonical ensemble where 
the coordinate partition function is given by eq.(|rT]). We thus used 



^ + — M(ri,...,r^) 



3Af/2-l 

e 



^+ ^M(fi,...,f^) 



as the statistical weight for the Metropolis algorithm. 

The MCE and CE Monte Carlo results coincide up to the statistical error for < < 
rjT, that is for 00 > ^ > — —0.19. In the MCE the gas does not clump at rj = rjc (point 
C in flg. 0) and the speciflc heat becomes negative between the points C and MC. In 
the MCE the gas does clump at ,^ ~ —0.52 , rj\.[(j ~ 1.33 (point MC in flg. |1|) increasing 
both its temperature and pressure discontinuously. We flnd from the Monte Carlo 
data that the temperature increases by a factor 2.4 whereas the pressure increases by a 
factor 3.6 when the gas clumps. The transition point rj'^jc Monte Carlo simulations 

is slightly to the right of the critical point rjMc predicted by mean fleld theory. The mean 
fleld yields for the sphere rjMc = 1-2598 . . .. 

In ref.[^ flnite corrections to the critical point rjMc are computed in mean fleld 
for the sphere. This flnite N corrections shift rjMc by +3.3% for = 2000. Since, rj1.jfj 
differs from rjMc by +5.6%, rjfjQ and rjMc are probably different critical points. 

As is clear, the domain between C and MC cannot be reached in the CE since cy > 
in the CE as shown by eg. (|5(]|) . 

We flnd an excellent agreement between the Monte Carlo and Mean Field (MP) results 
(both in the MCE and CE). (This happens although the geometry for the MC calculation 
is cubic while it is spherical for the MP). The points where the collapse phase transition 
occurs {rjT and rjj^jc) slowly increase with the number of particles A^. 

We verifled that the Monte Carlo results in the gaseous phase [rj < rj^) are cutoff 
independent for 10~^ > a > 10~^. 

As for the CE, the Gibbs free energy is discontinuous at the transition in the MCE. 
The transition is then of the zeroth order. We flnd from eq. (^) 



<l>(collapse) - ^{rjr) T^u , , n v ^ n 
— = jr- /(collapse) - /(r/^) ^ 0.7 > . 

gas gis 

where we used the numerical values from the Monte Carlo simulations. Notice that the 
Gibbs free energy increases at the MC transition whereas it decreases at the C transition 



[see eq.(|66|)]. 

Here again the two phases cannot coexist in equilibrium since their pressures and 
temperatures are different. 
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Microcanonical, gaseous phase, zi=-0.5 



one particle + 




Figure 3: Average particle distribution in the gaseous phase from Monte Carlo 
simulations with 2000 particles in the microcanonical ensemble for ^ = —0.5, t] = 
1.38, pV/[NT] = 0.277. 



We display in figs. |^-^ the average particle distribution from Monte Carlo simulations 
with 2000 particles in the microcanonical ensemble at both sides of the gravothermal 
catastrophe, i. e. rj = r]Mc- Fig. || corresponds to the gaseous phase and fig. § to the 
collapsed phase. The inhomogeneous particle distribution is clear in fig. ^ whereas fig. || 
shows a dense collapsed core surrounded by a halo of particles. 

The different nature of the collapse in the CE and in the MCE can be explained using 
the virial theorem [see eq. (|38D ] 

pV ^ U 
1 + 



NT NT 

When the gas collapses in the CE the particles get very close and U becomes large and 
negative while T is fixed. Therefore, may become large and negative as it does. 
We can write the virial theorem also as, 

1 1 
pV - -NT = - E . 
^2 3 



NT 



When the gas is near the point MC, E < is fixed and we have T > 0. Therefore, 
as well as U = E — 3 N T/2 cannot become large and negative as in the CE collapse. 
This prevents the distance between the particles to decrease. Actually, the Monte Carlo 
simulations show that < r > increases by 18% when the gas collapses in the MCE. 

Figs. ^ and |^ depict the average particle distribution from Monte Carlo simulations 
with 2000 particles in the canonical ensemble at both sides of the collapse critical point. 
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Microcanonical, collapsed, zi=-0.6 




Figure 4: Average particle distribution in the collapsed phase from Monte Carlo 
simulations with 2000 particles in the microcanonical ensemble for ^ = —0.6, t] = 

0. 43, pV/[NT] = 0.414. 

1. e. rj = rjc- Fig. |^ corresponds to the gaseous phase and fig. || to the collapsed phase. 
The inhomogeneous particle distribution is clear in fig. |^ whereas fig. |^ shows a dense 
collapsed core surrounded by a very little halo of particles. 

Notice that the collapsed phases are of different nature in the CE and MCE. The core 
is much tighter and the halo much smaller in the CE than in the MCE. 

Figs. 3 and 5 depict the average particle distribution for the gaseous phase in the MCE 
and the CE, respectively. In this phase, the MC simulations give identical descriptions 
for large N in both ensembles. [This important point will be further demostrated in sec. 
VI by functional integral methods] . The average configurations in figs. 3 and 5 describe a 
self-gravitating gas in thermal equilibrium within a cube. We may call it the isothermal 
cube by analogy with the well known isothermal sphere p0[1-||16||. 

6 Mean Field Approach 

Both in the microcanonical and the canonical ensembles the coordinate partition functions 
are given by 3A^-uple integrals [eqs. (|Tl|) and (|35D , respectively]. In the N ^ oo limit both 
3A^-uple integrals can be recasted as functional integrals over the continuous particle 
density as we see below. 
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Canonical, Gas Phase, eta = 1 .5 



one particle + 




Figure 5: Average particle distribution in the gaseous phase from Monte Carlo simulations 
in the canonical ensemble for 77 = 1.5 and = 2000 

6.1 The Canonical Ensemble 

We now recast the coordinate partition function e*^*^''-' in the canonical ensemble given 
by eg. (|35|) as a functional integral in the thermodynamic limit. 

A.^>1 f fjjpd^ ^-Nsa[pi.),a,v] (71) 



sc[pi.),d,r]] = -| y ^^^^^ P^^^ P^y^ + / Pi.^) logp(^) £xp{x) - ij 

where we used the coordinates x in the unit volume. The first term is the potential energy, 
the second term is the functional integration measure for this case (see appendix A). Here 
A^ p{x) stands for the density of particles. 

The integration over d enforces the number of particles to be exactly A^: 

d^x p{x) = 1 (72) 



That is, in the coordinates q (running from to L), the density of particles is 

A^ 
L3 



p{q) with J d^q — p{q) = N 
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Canonical, Collapsed Phase eta = 1 .53 



+ one particle 




Figure 6: Average particle distribution in the collapsed phase from Monte Carlo simu- 
lations with 2000 particles in the canonical ensemble for rj = 1.53, pV/[NT] = —14.44. 



6.2 The Microcanonical Ensemble 



Let us express the coordinate partition function in the microcanonical ensemble w{^,N) 
defined by eq. (pH]) in terms of the coordinate partition function in the canonical ensemble 
g^ivW defined by eq.(| 



In order to do that we use the Fourier expansion ||23 

r(A + i) du 



e{x) 



27r 



(73) 



We thus find from eqs.(0), {^) and that 



27r 



V 2 / A 27rz 



■log(Afr)) 



(74) 



where we introduced the integration variable rj = iuj/N and where 7 is an upward inte- 
gration contour parallel to the imaginary 77 axis. Using Stirling's approximation for the 
r function we find for ^ 1 up to irrelevant constants 



= L2m' 



■ log?? 
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Now, inserting the functional integral representation (fflD for the coordinate canonical 
partition function yields, 



J J 2 1X1 



(75) 



We thus find a functional integral representation in the micro canonical ensemble analo- 
gous to the canonical representation eq.(|7lD but with an extra integration (over 77) that 
constrains the value of the energy. 

The 'effective action' in the microcanonical ensemble takes thus the form, 



^ (76) 



6.3 The Grand Canonical Ensemble 

The partition function in the grand canonical ensemble can be written as 

00 

N=0 



(77) 



where z = stands for the fugacity and Z{N, T) is the partition function in the canonical 
ensemble given by eqs. (|28|) and (pT]). 

As shown in ref.[|], this grand canonical partition function can be expressed as a 
functional integral 



where 



Zgc{z,T) 



12 T 



^e.ff JV ^L2 



-eff 



An 



(7J 



(79) 



71 rj. 

Notice that the representation (|78D is exact while the functional integral representations 
in the microcanonical and canonical ensembles only apply for large number of particles. 
Rewriting eq.(^) in terms of the dimensionless variables (^ yields for the exponent 



Teff JV 



L 



d X 



2 ^<fix) 



where ^ = M L is of the order one (L^), since ~ z = ~ [see eq.(^)]. 

Since the exponent in the functional integral (|78D is proportional to L, the large volume 
limit is dominated by the stationary points (mean field approximation) 



^0) 



We expand around the saddle point $s(a^) changing to a new functional integration Y{x) 
variable as follows, 

$(f) = $^(f) + Y{x) . (81) 
Keeping in eq. ([78|) quadratic terms in Y{.) yields, 



ZGciz,T) = e 



1 + 

(82) 
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where the Gaussian integral over Y{.) gives a factor of order L° [see paper II]. 

We recall that the saddle point method applies while all eigenvalues of the quadratic 



form in the exponent of eq.(|82D are positive. Therefore, the determinant of the quadratic 
fluctuations is positive. The determinant vanishing or changing sign indicates the presence 
of zero or negative eigenvalues. In such a case the system is no more on a stable phase 
but on a metastable or unstable phase. The free energy gets an imaginary part in such 
metastable or unstable situations. 

The average number of particles in the grand canonical ensemble is given by 

1 ^ .r dlogZcc 



N 



Zgc n=q 



9 log z 



We thus find in the mean field approximation, 



N 



eff 



(fx e*=(^) 



Therefore, using this and eq.([75|) we can express in terms of rj where we denote N 
as N to avoid cluttering of notation. 



and the fugacity results 



N I 2ti \3/2 1 



(83) 



^4) 



We again see that z L ^ in the GCE. 

Integrating eq.(PD|) over the unit volume yields 



j V$s(f) ■ ds = -Anr] 



15) 



where we used eq. 

We find for the free energy P], 



NT 



F = -T log Zgc + NT log z = Fo+ ^-^K{r]) - NT log + 0{N^) , (86) 



where we used the grand canonical partition function (|8^ ) evaluated at the stationary 
point. 



log Zgc = N 

and z is given by eq. (|8^) with 

Jod^x <l>s{x) e*=(^) 



^7) 



Fq is given by eq. 



C{v) 



and CM = d^xe"^'^^^ . 
Jo 
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We easily calculate the mean value of the potential energy in the mean field approxi- 
mation 

<U>=-T = -^KM (89) 

(7 log G 2 

Combining the two expressions for the entropy 

E-F fdF\ 
S = _ aud 5=-(^)^. (90) 

yields 

S = So-N[Kiv)-\ogCir])] (91) 
and the first order differential equation 

T] K'{ri) + K{ri) = 27]^ log C{7]) . (92) 

drj 

The boundary conditions K{Q) = 0, C(0) = 1 ensure the ideal gas limit r/ = 0. 
The pressure takes the form, 



p^_(dF\ NT 



+ 0(Ar°) . (93) 



These equations guarantee in addition that the virial theorem (^) holds. 

6.4 Saddle point evaluation in the canonical ensemble 

The functional integral in eq.([7T|) is dominated for large by the extrema of the 'effective 
action' sc[p(.), a, 77], that is, the solutions of the stationary point equation 

, /^x f d^y Ps{y) 

\ogps{x) - T] / — — — = as , (94) 
J \x — y\ 

a = id is a Lagrange multiplier enforcing the constraint (|72D . 
Applying the Laplacian and setting (f){x) = logps(a;) yields, 

V^0(x) + Attt] e'^(^) = , (95) 

This equation is scale covariant |Q. That is, if 0(x) is a solution of eq. (|95|) , then 

(j)x{x) = (t){Xx)+\ogX^ (96) 



where A is an arbitrary constant is also a solution of eq . (pSf) . For spherically symmetric 
solutions this property can be found in ref.|p]. 

Integrating eq. (|95D over the unit volume and using the constraint (^) yields 

V(f){x) ■ ds = -Anri (97) 

where the surface integral is over the boundary of the unit volume. 

Comparing eqs.(^)-(^) with ( P^D and ( P7| ) shows that the grand canonical and canon- 
ical stationary points are related by 

$,(f) = 0(f) + log C(r7) . (98) 
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Eq.(p2D can then be written as 

^ , rr-. 4^ |/'d3^^[i<^V20+47r,?e^(S)]-27r,?logC(r,)| 

^Gc[z, 1 ) — e ^ J 



X 



(99) 



where we used eqs.(|32|), (0), (||) , (|93) and (H). 



We have taken the zero cutoff hmit in eqs.(p^)-(p5|). The mean field equations turn 
to be finite with regular solutions in such limit. This can be understood from our 
perturbative calculation in sec. III. A. All potentially divergent contributions at zero 
cutoff are suppressed by factors and therefore disappear in the N = oo limit. Hence 
one can set the cutoff to zero in the mean field approximation. 



In order to evaluate the functional integral in eq.(|7TD by the saddle point method we 
change the functional integration variable as follows, 



p(x) = ps{x) + Y{x) 



a = as + Vo 



where Ps{x) and obey eq.(0). We can expand the exponent to second order as 

sc[p{.), a, ri] - sc[ps{.), a„ r]] = [¥{.), yo] + O (r^ y^) 
where we use that 



(100) 



:ioi) 



6sc 



6p{x) 



ds 



c 



da 



p=ps,a=ag 



and 



sg)[r(.),l/o] = ^ I d'xd'yY{x)Yiy) 



Notice that 



5p{x)5p{y) 
do? 



p=ps,a=as 



+ yo (Tx Y{x) 



Sp{x)da 



p=Ps,a=as 







We evaluate exphcitly the second derivatives from eq.(^) with the result, 

6^sc S{x-y) 7] , 



6p{x) 6p{y) p{x) 



\x - y\ 



Sp{x) da 



Therefore 

.(2) 



1 

2 



d'^x 



Y'^{x) 1] f d^x d^y 
p{x) 2 J \x — y\ 



Y{x) Y{y) - yo / d^x Y{x) (102) 



Inserting eqs.( |100| ) and ( 101|) into eq.(|7ll) yields 

J J DYdyo e-^4''in),J/o] 



l + O 



(103) 



where s{r]) = sc[ps(.), fls, stands for the value of the exponent at the saddle point. 
Terms of order higher than quadratic in sc[p{-),a,ri] contribute to the 1/A^ corrections. 
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The Gaussian functional integral (|103| ) can be exactly computed in terms of the func- 
tional determinant of the quadratic form ( |102| ) [see paper II]. It gives a result of order one 
(iVO). 

In the mean field approximation we only keep the dominant order for large A^. There- 
fore, only the exponent at the saddle point accounts and according to eq.(|3l[) we find for 
the free energy 

F = Fq + NTs{7]) + 0{N^) 

Hence, in the mean field approximation, the function f{r]) is given by 

fMFiv) = 1 + , (105) 
6 drj 

From eq.([7TD we can compute s{ri) in terms of the saddle point solution as follows 

siv) = sc[ps{-),as,rj\ = [ ^.^ ^ ^1. ps{x) ps{y) + / d^x ps{x) logp,(x) . (106) 

1 J \x — y\ J 

Using eq.(^4|) we find an equivalent expression that will be useful in paper II, 

s(ry) = ^ + 1 1 0(f) e'^(^) d\ . (107) 

6.5 Saddle point evaluation in the microcanonical ensemble 



The extrema of the 'effective action' ([76|) dominate the microcanonical partition function 
([75|) in the large limit. Extremizing eq.(|7B|) with respect to p(.) and d gives again 
eqs.(0) and (|72|), respectively. 

An extra equation follows by extremizing the 'effective action' ([7B|) with respect to rj: 



3 1 /■ d?x d^y 
2'qs 2 J \x-y\ 

Going back to dimensionful variables this equation takes the familiar form 

3 Gm'^ rd^qd^q' 



2 2 J \a — Q \ 



q-q' 

That is, eq. (|108|) enforces the fixed value of the energy in the microcanonical ensemble. 

Therefore, the stationary point equations in the microcanonical and canonical ensem- 
bles are identical. Both ensembles yield the same results in the N ^ oo limit in their 
common region of validity. We derive the domain of validity of the mean field approach 
for each of the three statistical ensembles in paper II. That is, the regions where all fluc- 
tuations around it decrease its statistical weight within their common region of validity. 

In order to evaluate the functional integral for the microcanonical partition function 

m 

J J 2 TTIi 

30 



we expand the 'effective action' SMc[p{.-)iO.^'n\ around the stationary point Ps{-),CLs,r]s to 
second order. This gives 



(110) 



where Y{.) and yo are defined by eq. (|100|) and we set rj = rjs + rj. The second order piece 



of the 'effective action' takes now the form 

.(2) r\^M -1- (2)rvM 1 - f d^x d?y 3 



stlc[y{-),y^.fl] = s'^\Y{.),y,]-f, / — — ^ p,(f) Y{y) - —-f,' . (Ill) 

J \x y\ 

The Gaussian functional integral in eg. ( |110| ) yields a contribution of order one {N^) [see 
paper II]. The dominant (mean field) contribution, e~^^^'^\ exactly coincides with 
the mean field result in the canonical ensemble [eq.( p.03| )] Therefore, the canonical and 
microcanonical ensembles yields identical physical magnitudes and the same equation of 
state in the mean field limit. 

6.6 Spherically symmetric case 

We shall consider here the spherically symmetric case where eq.(^) takes the form 

where we work on an unit volume sphere instead of an unit volume cube as in eq.(^. 
Therefore, the radial variable runs in the interval 

It is more convenient to introduce a new radial variable 



r = R 



3 



such that < r < 1. 

The saddle point equation ( |112| ) takes then the form 



^ + 2^ + 4vr,-e^(^) = 0. (113) 
ar^ r ar 

where 

/47r\i/3 3 

rj^ = n ( — ) = 1.61199... and e-^^") = e'^(^) — . (114) 
V 3 / 47r 

In order to have a regular solution at r = one has to impose 

0'(O) = . (115) 



Otherwise, the second term in eq. (|112|) diverges for r ^ 0. 



In the spherically symmetric case, the constraint (|9y|) becomes 

0'(1) = -r/^ • (116) 
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Figure 7: t]^ as a function of the uniformizing scale variable log A according to eq. (|120|) . 
Notice the maximum of rj^ as r^^ = 2.517551 . . .. The region beyond the point MC, 
(In Xmc = 3.53698 . . .), is unphysical as we discuss in paper II 
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Using the scale covariance ( ^6]) we can express (f){r) as 

/ ^2 \ 



log 



4 711]^ 



where 



X"(A) + ^ X'(A) + e^W = , x'(0)=0 



:il7) 



(118) 



This equation is invariant under the transformation: 

A^Ae" x(A)^x(A)-2a, (119) 

where a is a real number. Hence, we can set x(0) = without loosing generality. 
x{x) is independent of r]^, and A is related to r]^ through eq. (|116|) 

A X'(A) = -r/« . (120) 

Since A and r]^ are always positive, x(A) is a monotonically decreasing function of A. 
Eq. (|118|) can be easily solved for small arguments as 

2 4 

Hence, in the dilute limit eq. (|120|) relating r]^ with A gives 



V 



R 



X2 \4 

--- + 0(A«). 
3 30 ^ ' 



[121) 



For large argument, the solution of eq.( |118D takes the asymptotic form| 

'V7 



X{x) = log — + — cos 



X 



■ log X + B 



l + O 



:i22] 



where A and B are numerical constants. Using eq. (|120|) this gives for rj 



R 



(123) 



where C and D are constants related to A and B. By numerically solving eq.( |118| ) we 
find 

C = 1.667... . 

It must be noticed, however, that the mean field solution is unphysical for A > Xmc = 
34.36361 ... as we shall see in paper II. Anyway, we see from fig. ^ that r/^ approaches 
very fast its asymptotic behaviour ( |123|) for log A > 2. 
We plot in fig. |^ xi^iv^)) a function of rj^. 

In the spherically symmetric case the integral over the angles in eq . (|9^) is immediate 
with the result, 



tts + 4:711] 



1 

r Jo 



r'^ dr' e^('^') + / r' dr' e'*'^''^ 



(124) 
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Deriving with respect to r yields, 



dr 



Jo 



This again shows that 0(r) is a monotonically decreasing function of r [see above, eg. ( |120| )]. 
Setting r = 1 in eq.( |124| ) leads to the relation 



'I) = a, 



+ Attt]^ dr e'^^'^ 
Jo 



Using now the constraint (^) allows us to compute the Lagrange multiplier a at the 
saddle point 

(125) 



as = 0(1) - 1] 
The particle density in MF is given by 



R 



p{r) = e 



^(^"^ , < r < 1 . 



Since xW monotonically decreases with A, the particle density monotonically decreases 
with r for fixed rj^. 

Let us now compute s{ri^) [the exponent in eq.(^) at the saddle point] for the spher- 
ically symmetric case. We find from eq.( |107|) 



s{ri 



1 

2 

log 



1) - +2n I r'^ dr (t){r) e 



Airrj^ 



2 2 A r/^ Jo 



x'^ dx [x'(x)]^ 



(126) 



where we integrated by parts and used eqs. (|117|) -( |120| ). 

The integral in the r.h.s. of eq.( |126| ) can be computed in closed form [see appendix B] 
with the result, 



s{ri^) = log 



+ ^(A) + 3-r^«-4e^(^) 



Airrj^ 



n 



R 



Inserting now s{ri^) into eq.( |105|) and using eqs. (|118|) -( |120|) yields after calculation 



s(r/^) 



A2 



3[l-/MF(r7^)]-r7^ + log 



(127) 



3/a/f(?7^ 



Notice that fMpiv^) well as the other physical quantities are invariant under the 
transformation (11191) as it must be. 



It follows from eqs.( |118| ), ( |120| ) and ( |127| ) that fMpiv^) obeys the first order non-linear 
differential equation 



r7^(3/A/F - l)fMFiv'') + (3/a/f - 3 + r]'')fMF = . 



:i28l 
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' pV/[NT] = f(eta)'vs. eta 
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Figure 9: fuFij]^) = PV/[NT] as a function of r]^ in the MF approximation [eq.( |128| )]. 
fhiFijI^) has a square root branch point at rj^. The points GC, C and MC indicate 
the transition to the collapsed phase for each ensemble (grand canonical, canonical and 
microcanonical, respectively): rjQ(j = 0.797375 . . . , 77^ = 2.517551 . . . , tj^q = 2.03085 . . . 
(notice that rj^jQ is in the second Riemann sheet). Since E/[3NT] = fMpij]^) — |, this 
plot also shows the energy per particle as a function of r]^. Furthermore, the particle 
density at the surface is proportional to fMpijI^) [eq. (|129|) ]. 
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which reduces to an Abel equation of first kind[p2 

We thus find that in the mean field approximation all thermodynamic quantities follow 
from the resolution of the single first order non-linear differential eq. ( |128| ) with the initial 
condition /mf(0) = 1. 

Integrating eq.( |128| ) with respect to 77^ yields, 

~ ^'"^^"^^^ = 3[/,,^(r^^) - 1] + r/^ - log fMPiv"") 
Further useful relations follow from eqs.( |117| ) and ( |127| ) 



0(1) = log 



3/AfF(r7^) 



An 



P(l) = ^ fuFiv'') ■ 
4 vr 



(129) 



That is, the particle density at the surface (r = 1) is proportional to fuFij]^)- 

We can then write the different physical magnitudes in the MF approximation as 



pV 

Irf 

F-Fq 

NT 
S-Sq 
N 
E 

iVT 



3[l-/A/F(r/^)]-^'' + log/MF(^^) 
6[/MF(r7^)-l]+r7^-log/AfF(r/^) 



(130) 



where we used eqs.(|41| 



1), (H) and 

We derive in appendix C the properties of the function fMpiv^) from the differential 
equation (|128|) . One easily obtains for small rj^ (dilute regime). 



fMpiV 



R\ 



R\2 



175 



These terms exactly coincide with the perturbative calculation in the dilute regime for 
spherical symmetry [see eq. (P^D , ( pTj ) and ( |114D ]. 

We plot in fig. 1 fMpiv^) ^ function of rj^ obtained by solving eq. (|128| ) by the 
Runge-Kutta method. We see that fMpiv^) is a monotonically decreasing function of 
T]^ for < T]^ < 1]^. At the point r]^ = t]^, the derivative fMpiv^) takes the value —00. 
It then follows from eq. (|128|) that 

f'MpiVc) = I ■ 



At the point the series expansion for fupij]^) in powers of rj^ diverges. Both, from 
the ratio test on its coefficients and from the Runge-Kutta solution, we find that 



f]^ = 2.517551 



:i3i) 



From eq.( |128| ) we find that fMp{v^) ~ \ has a square root behaviour around rj^ = rj^: 



JMP[V ) - 



2{^§ - 2) 



1]^ + 



2 - 1) 



7Vc 
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Figure 10: The entropy per particle minus the ideal gas value as a function of rj^ in the 
MF approximation [eq.( |130| )]. 



Inserting the numerical value ( |131| ) for rj^ yields, 



fMFiv'') ' I + 0.213738 . . . - r/« + 0.172225 . . . (ry^ - r;^) + O [(r^^ - v^^'' 

(132) 

We see that fupiv^) becomes complex for t]^ > t]^. Recall that in the Monte Carlo 
simulations the gas phase collapses at the point r]^ < rj^. 



From eqs.(^), we plot pV/[NT], S/N and as a function of r]^ in figs, g |TU 

and ITT], respectively. 



The points GC, C and MC correspond to the collapse phase transition in the grand 
canonical, canonical and microcanonical ensembles, respectively. Their positions are de- 
termined by the breakdown of the mean field approximation through the analysis of the 
small fluctuations [see paper II]. 

fMpijI^) is a multivalued function of rj^ as well as all physical magnitudes [see eq. (|130|) ]. 
As noticed before, the CE only describes the region between the ideal gas point, rj^ = 
and C in fig. 1. The MCE goes beyond the point C (till the point MC) with the physical 
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magnitudes described by the second sheet of the square root in eqs. (|132| ) (minus sign) 
We have near C between C and MC, 



= - - 0.213738. . . 



77^ + 0.172225... [t]^ 



i?.\3/2" 



The function fMpiv^) takes its absolute minimum at r]^ = rj^i^ — 2.20731 ... in the 
second sheet where fMpiVmin) = 0.264230 . . .. 

Since fuFij]^) < \ imphes that the total energy is negative [see eq. ( |130| )], the gas is 
in a 'hounded state' for rj^ beyond = 2.18348 ... in the first sheet. 

Since xW and ri{X) are single- valued functions of A, /a/f(^^(A)) defined by eq.( p.27| ) is 
also a single- valued function of A. That is, A is the uniformization variable. All physical 
magnitudes are single- valued functions of A. On the other hand, A is an infinite- valued 
function of 7]^ as one sees from fig. |^ and eq.( p.23| ). That is, fMpiv^) has an infinite 
number of Riemann sheets. However, only the first two sheets are physically realized. 
The rest are unphysical. A plot of fuFij]^) including all sheets produces a nice spiral]^ 
converging towards t]^ = 2 , fMpiv^) = 1/3 for A = oo as follows from eqs. (|122|) , (|123|) 
and (11271) . 



A induces a scale transformation in coordinate space as we see in eq. (|1 1 7| ) whereas 



f]^ plays the coupling constant [Recall that rj^ is proportional to Newton's gravitational 
constant]. 

The variation of r]^ with respect to A yields the renormalization group equation 

A^ = r/^[3/MH^^)-l] 

where we used eqs. (|118|) , ( |120|) and (|127|) . Here r]^ fuFijI^) — 1] plays the role of the 
renormalization group beta function. We see that it has two fixed points at r/^ = and 
at f]^ = TjQ. [See fig. ^ where the running of r]^ with A is exhibited]. 
We find from eqs. (|121|) and ( |132|) near these fixed points 



R A-^o 

^ = y 

^R A^A. ^« _ 2) _ y^^y 

c 

where the coefficient has the numerical value -^^^-ft — - = 0.0085515 . . .. 

6.7 Canonical vs. Grand Canonical Ensembles in the Mean 
Field Approximation 

We have seen that the stationary point equations and their respective solutions are closely 
related in the canonical and grand canonical ensembles [eqs.(^0|)-(^) and (|97|)-(|98D]. 

Let us now show that physical quantities obtained from both ensembles do coincide 
in the mean field approximation. 

From eq.(|88|) and (|98D we find that 

K{r]) = J (f){x) e'^(^) d^x + logC(r/) . (133) 
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Figure 12: {cv)mf as a function of r]^ from mean field eq. (|137|) . Notice that {cv)mf 
diverges at the point C, that is for rj^ = 2.517551 .. . 



[Recall that / e'^(^) d^x = 1]. 

In the spherically symmetric case this integral takes the form 

2 

= 6[l-/MF(r7^)]-r7^+log 



47r dr 0(r) e'^^''^ = 0(1)+^? dr 

Jo 7]^ Jo 



dr 



R\ 



where we integrated by parts and used eqs.( |117[ ), ( |120D and appendix B. 
From eqs.( |133| ) and ( |134|) we find 

3/m£(?7 
47r 



(134) 



Kir]^) - logC(r/«) = 6 [1 - fMpiv"")] - v"" + log 



R\ 



R\ 



Inserting this result into the linear differential equations (0) leads to the solution, 

Civ') = ^f^, and Kiv^) = Q[l-fMFiv'')] 
3 fMF{ri^) 

We then find from eqs.(||), (|12|) and ( [129D that 

logC(r7) = -as . 
Combining eq.(|l3|) with cqs.®, (ID-dlD and (| 



(135) 



(136) 

shows that the canonical and the 
grand canonical ensembles yields identical physical magnitudes (pressure, energy, en- 
tropy, free energy, specific heats, compressibilities, speed of sound) and the same equation 
of state in the mean field approximation. 
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Figure 13: {ht)mf as a function of r] from mean field eq.( |139| ). Notice that {kt)mf 
diverges at t]^ = rj^ = 2.43450 . . . 

The thermodynamical potential]^, 

Q = -T log Zgc = N[3 fMEiv"") - 2] 
is not equal to —PV. That is, here Q ^ —PV and we have instead 

Q + PV = 2NT[l-fMF{v'')] 



This relation is analogous to eg . (^81) . differs here from —PV since for the self-gravitating 
gas we have N L instead of the usual relation N L^. 



7 Specific Heats, Speed of Sound and Compressibil- 
ity 



The specific heat at constant volume in the mean field approximation takes the form 

(137) 



7 — 2 

{cv)mf = 6 fMFiv^) -7; + V^ + T7 — rifT 

2 3 fMF[V ) 



where we used eqs. (|i9|) and (|128|) . 

We plot in Fig. eq. (|137D for {cv)mf as a function of rj. We see that {cv)mf 
increases with i] till it tends to +oo for t]^ f rj^. It has a square-root branch point at the 
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Figure 14: cp as a function of r] from mean field eq.( p.43| ). Notice that cp diverges at 



j^R = = 2.43450 . . . 



point C. In the stretch C-MC (only physically realized in the microcanonical ensemble), 
{cv)mf becomes negative. We shall not discuss here the peculiar properties of systems 



with negative Cy as they can be find in refs.jll, 12, 16 



From eqs. (|132|) and ( 137 ) we obtain the following behaviour near the point C in the 
positive (first) branch 



(cv)mf = 0.80714 ...{vc- v"")'^'^ " 0.19924 . . . + 0{^r]^ - r]^ 
and between C and MC in the negative (second) branch 

(cv)mf ^ -0.80714 ...(vc- v"')''^^ " 0.19924 ... + 0{ 
Finally, {cv)mf vanishes at the point MC rj^^fj = 2.03085 . . .. 



Vc - V^) 



The isothermal compressibility in mean field follows from eqs. (p^ and ( p.28| ) 



[kt)mf 



VuFm 



1 + 



QfMFiv^) - 



(138) 



(139) 



We plot {ht)mf in fig. 0. We see that {kt)mf is positive for < r]^ < tiq = 2.43450 . . . 
where {kt)mf diverges. The point rj^ is defined by the equation 



QfMF{Vo)-Vo=0. 



(140) 
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Figure 15: The speed of sound squared at the surface divided by the temperature, vl/T, 
as a function of rj from mean field [eq.( |143|) ]. Notice that v1/T takes the value 11/18 at 
the critical point 1] = rjc and becomes negative beyond rj^ = 2.14675. . . in the second 
sheet. 



We find from eqs. (|12|) and ( |140| ) that 



fMpiv^) = -\ (141) 



{kt)mf diverges for r]^ ~ rj^ as 



{i^t)mf is negative for < r]^ < r]^ and exactly vanishes at the point C. {ht)mf 
then becomes positive in the stretch between C and MC only physically realized in the 
microcanonical ensemble. 

Notice that the singularity of {k.t)mf at rj^ = rj^ = 2.43450 ... is before but near the 
point C. It appears as a preliminary signal of the phase transition at C. 7]q is probably 
the transition point tjt seen with the Monte Carlo simulations (see fig. 1). (Recall that 
rjT ~ 1.515 corresponds to r^^ ~ 2.44). 

It is easy to understand the meaning of a large compressibility. From the definition 

(0) 

f^.KrSp=-.r^. (142) 
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A large compressibility implies that a small increase in the pressure {5p <^ NT/V) pro- 
duces a large change in the density of the gas. That means a very soft fluid. 

For negative compressibility, eq. (|142|) tells us that the gas increases its volume when 



the external pressure on it increases. This is clearly an unusual behaviour that leads to 
instabilities as we shall see below. 

The specific heat at constant pressure in the mean field approximation takes the form 

3 24 (t]^ - 2) fMFiv^) 



where we used eqs.(p3D and ( |128| ). We plot {cp)mf in fig- 0- We see that (cp)mf 



is positive and grows with r]^ till it diverges at the same point where {ht)mf diverges 
ri^ = n^ = 2.43450 .... That is, 

— V Vo ~ V 

{,cp)mf becomes negative for tjq < r]^ < r]^. It keeps negative in the C-MC section till 
the point t]^ = 7]^ = 2.14675 . . . where it becomes positive. The point rji is defined by 
the equation 

2 

The speed of sound squared at the surface in the mean field approximation takes the 
form 



24 /^^(r^f) + (4r^f - 19) /^^(r/f ) + ^ = . (144) 



T 3 



^1 3/mH^^) + V-2 

QpMFiv'') + {v''-'i)fMFm + i 



(145) 



where we used eqs.(|^ and ( |128| ). We plot ^ as a function of in fig. |T3|. We see that 
^[t]^) is positive and decreasing with rj^ in the whole interval between r]^ = and C. 

2 2 

At the point C it takes the value if^{ric) = 11/18. Then, ifr{r]^) decreases between C and 
MC becoming negative at t]^ = 2.14675... in the second sheet where it vanishes. Notice 

2 

that jriv^) {cp)mf vanish at the same point rj^ defined by eq.( p.44| ). 

Vg<0 indicates an instability. That is, small density fluctuations grow exponentially 
in time instead of propagating harmonically. It is remarkable that Vg becomes negative 
at rji = 2.14675... in the second sheet before but near the MC critical point rjfjc = 
2.03085... in the second sheet. Somehow, the change of sign in v"^ announces the MC 
critical point. 

ifiv^) tends to —oo for rj^ | Vmc- Notice that the denominator in eq. ( |145| ) exactly 
vanishes at rj^ = rjf^Q [see Table 1]. 

The adiabatic compressibility ^5 is not here an independent quantity. We find from 
eqs.(ll), O, O and (|143D, 

cv 3 12/|,^(r7«) + (2r7^-ll)/MF(r/^) + l 



Ks = Kt 



CP fMFm 48 fl^^r^R) + (8 r^^ - 38) /mf^ + 



That is, 
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POINT 


A 




Defining Equation 


Imf^ 


PHYSICAL MEANING 


GC 


1.7772... 


0.797375 . . . 


2 - 3 77gc fMFiVGc) = 


0.836076 . . . 


Collapse in the GCE. 


3 


3.38626 . . . 


1.73745... 


3-r?^ + x(A) = 


0.622424. . . 


Energy density 
vanishes at r = 0. 


2 


4.73739 . . . 


2.18348... 


2 fMF{v2) -1 = 


0.5 


Total Energy vanishes. 





6.45077... 


2.43450 . . . 




0.40575 . . . 


kt and cp diverge. 


C 


8.993195... 


2.517551... 


1 - 3 fMFivE) = 


1/3 


Collapse in the CE. 
cv diverges. 


Min 


22.5442 . . . 


2.20731 . . . 


fMpiVmin) = 


0.264230. . . 


Minimum of 

pV/[NT] 

iu the gas i)hase 


1 


25.7991 . . . 


2.14675 . . . 


Wmf{v^) - (38 - 8r/f ) /mf«) + 
+77f = 


0.265290... 


Vg and Cp vanish. 


MC 


34.36361 . . . 


2.03085 . . . 


Wmf{v^c) - (11 - 277^c) /mf«c)+ 
+1 = 


0.273512... 


Collapse in the MCE. 
Cv vanishes. 



8 Discussion 



We have presented here a set of new resuhs for the self-gravitating thermal gas obtained 
by Monte Carlo and analytic methods. They provide a complete picture for the thermal 
self-gravitating gas. 

Contrary to the usual hydrostatic treatments 0, we do not assume here an 
equation of state but we obtain the equation of state from the partition function [see 
eq.(ET[)]. We find at the same time that the relevant variable is here 7]^ = Gm^N/lV^^^T]. 
The relevance of the ratio Gm'^ /[V^^^T] has been noticed on dimensionality grounds 
However, dimensionality arguments alone cannot single out the crucial factor in the 
variable t]^. 

The crucial point is that the thermodynamic limit exist if we let ^ oo and ^ oo 
keeping r]^ fixed. Notice that r] contains the ratio N V~^^^ and not N/V. This means 
that in this thermodynamic limit V grows as A^'^ and thus the volume density p = N/V 
decreases as ~ N~'^. rj is to be kept fixed for a thermodynamic limit to exist in the same 
way as the temperature. pV, the energy E, the free energy, the entropy are functions of 
1] and T times N. The chemical potential, specific heat, etc. are just functions of t] and 
T. 

We find collapse phase transitions both in the canonical and in the microcanonical 
ensembles. They take place at different values of the thermodynamic variables and are of 
different nature. In the CE the pressure becomes large and negative in the collapsed phase. 
The phase transition in the MCE is sometimes called 'gravothermal catastrophe'. We find 
that the temperature and pressure increase discontinuously at the MCE transition. Both 
are zeroth order phase transitions (the Gibbs free energy is discontinuous). The two 
phases cannot coexist in equilibrium since the pressure has different values at each phase. 

The parameter rj^ [introduced in eq.(^)] can be related to the Jeans length of the 
system 

dj = \^^^ , (146) 
V m y/G m p 

where p = N/V stands for the number volume density. Combining eqs. (|32D and ( |146|) 
yields 

We see that the phase transition in the canonical ensemble takes place for dj ~ L. [The 
precise numerical value of the proportionality coefficient depends on the geometry]. For 
dj > L we find the gaseous phase and for dj < L the system condenses as expected. 
Hence, the collapse phase transition in the canonical ensemble is related to the Jeans 
instability. 

The latent heat of the transition (g) is negative in the CE transition indicating 
that the gas releases heat when it collapses [see eq.(^)]. The MCE transition exhibits 
an opposite behaviour. The Gibbs free energy increases at the MCE collapse phase 
transition (point MC in fig.|lD whereas it decreases at the CE transition [point T in fig. 
|T], see eq.(p^]. Also, the average distance between particles increases at the MCE phase 
transition whereas it decreases dramatically in the CE phase transition. These differences 
are related to the MCE constraint keeping the energy fixed whereas in the CE the system 
exchanges energy with an external heat bath keeping fixed its temperature. The constant 
energy constraint in the MCE keeps the gas stable in a wider domain and makes the 
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collapse transition softer than in the CE. Notice that the core is much tighter and the 
halo much smaller in the CE than in the MCE [see figs. ^ and |^. 
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A Functional integration Measure in the Mean Field 
Approach 



We follow the derivation of ref. for the functional integral measure. We want to recast 



JO 



^ri u{ri,...,fM) 



(147) 



1=1 



as a functional integral in the large limit. 

We start by dividing the domain of integration (of unit volume) into M cells. Each 
cell is of volume Ur and contains kr particles with 1 < r < M. Therefore, 



M 



M 



J2kr = N , ^ = 1 . 

r=l r=l 

We can thus rewrite the multiple integral (|147|) as follows: 

/ M \ ATI M 

ki,...,kM V r=l I \.\.r=\'^r- r=\ 



where pT I 



1 11 

J = ^ k^i Vry + - ^ kr K-,rH ^ ^r' ^r" [< K-,r' K-,r" > — < Vr,r' > < K-,r" >] + • 



and 



Vr r' 



1 7] 



1 fl 



UJr UJr' N Jo Jo Ti — r2 



Assuming \/N <^ujr < N^"^^^ one can neglect in J terms quadratic and higher in K,r'|^ 
The particle density is defined as 

M u 

Therefore, we can write the sums over r as integrals in the following way 



— kr kr' Vr^ 



1 rl 



2N Jo Jo n — r2 
Using Stirling's' formula one finds that 



P(ri) p(f2) • 



M ( \kr 

n[UJr) N- 
k I 



M 

n 



1 
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-N J (fix p{x) log[p{x)/e] 



Collecting all terms yields, 

M f \k 3 3 

r=l ^f' 

whereas the constraint in the number of particles takes the form 
and finally, 

^^>i Jdp e^/Ti#''(^^)''(^^-^/'^^^^(^^ 5 (/ d'xpix) - l) 

Replacing the Dirac delta by its Fourier representation 

^8 [j d\p{x) - = I ^ ^.7Va(/d3.p(.^)-l) 



yields eq.([7TD. 

B Calculation of the saddle point 

We prove in this Appendix that the integral 

/(A) = dx [x'{x)]^ (148) 



takes the value 

/(A) = Xr^R (Q-r]^) -2 A^ e^^^^ (149) 
Here x{^) is a regular solution of eg .( 11181) in the interval < x < A fulfilling the relation 

We start by computing the derivative of /(A) in two ways. According to the definition 

(H) 



dX 



Then, we compute the derivative of eq.( |149| ) with respect to A and use eqs. ( |1 18| ) and 



11201). We find after calculation that both results coincide. 



Finally, we observe that both eqs.( p.48| ) and ( |149| ) vanish at A = 0. Therefore, eq. ( |149| ) 
is valid. 



C Abel's equation of first kind for the equation of 
state 

In the mean field approximation the equation of state for spherical symmetry satisfies the 
first order differential equation ( |128|) 

r/^(3/MF - l)fMFiv'') + i^fMF - 3 + r/^)/MF = . (150) 
with the boundary condition fupiS^) = 1- 
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We can solve eq.( |150| ) in power series in rj^ around the origin, 

oo 

fuFiv) = 1 + E /n ^" 



(151) 



n=l 



Inserting eq.(151) into eq.(150) yields the quadratic recurrence relation 

k=2 



fn 



2n + 3 



for n > 2 . 



where /i = -i. 

We find from this recurrence relation, 



/2 



175 



1575 



991 



3031875 



All coefficients fn are negative rational numbers for n > 1. They decrease very fast with 
n as 



fn 



n>l 



0.0956678 , 



l + 0{- 

n 



This formula reproduces the large orders of the expansion of \JriQ — rj^ describing the 
behaviour of fpiriv) ^^ar r/g [see eq.( |132[ ) and ref.[Q] 



Vc ~ 



1 /^gr(n-i) f^^y 



2 V vr — n\ 



Notice that 



1 /^r(n-i) „>i 0.447594... 



i + o{- 

n 



2 V vr n\ 

and that 0.213738 . . . x 0.447594 . . . = 0.0956678 . . .. 

The power series ( p.51[ ) thus has a radius of convergence 77^ = 2.517551 .. .. The 
singularity of fupij]) nearest to the origin is thus the critical point. 
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